#include <stdio.h>
#include "lattice.h"
#include <math.h>

void readdata( float T23[NFILES][NCOOLINGSTEPS][NT], 
			float T24[NFILES][NCOOLINGSTEPS][NT], 
			float T34[NFILES][NCOOLINGSTEPS][NT]);
float average(float *arr, int L);
float error  (float *arr, int L);
void makeOutput(float Correlator[NCOOLINGSTEPS][NT], float Error[NCOOLINGSTEPS][NT]);

//===============================================================================================
int main(int argc, char *argv[])
{
	float T23[NFILES][NCOOLINGSTEPS][NT], T24[NFILES][NCOOLINGSTEPS][NT], T34[NFILES][NCOOLINGSTEPS][NT];
	float Correlator[NCOOLINGSTEPS][NT];
	float Error[NCOOLINGSTEPS][NT];
	
	readdata(T23, T24, T34);

	int i,dt;
	for(i = 0; i < NCOOLINGSTEPS; i++)
	for(dt = 0; dt < NT; dt++)
	{	
		float tempcorr[NFILES];
		for(int k  = 0; k < NFILES; k++)
		{
			float avsq=0, av=0;
			for(int j1 = 0; j1 < NT; j1++)
			{
				int j2 = (j1+dt)%NT;
				avsq += T24[k][i][j1]*T24[k][i][j2];
				av += T24[k][i][j1]; 
			}	
			tempcorr[k] = avsq/NT - av*av/NT/NT;
		}
		Correlator[i][dt] = average(tempcorr, NFILES);
		Error[i][dt] = 1/sqrt(NFILES-1)*error(tempcorr, NFILES); 
	}

	makeOutput(Correlator, Error);
}
//===============================================================================================

void readdata(float T23[NFILES][NCOOLINGSTEPS][NT], 
			float T24[NFILES][NCOOLINGSTEPS][NT], 
			float T34[NFILES][NCOOLINGSTEPS][NT])
{
	FILE *input;
	input = fopen("out01", "r+b");
	if(input==NULL)
	{
		cout << "error reading file\n";
		exit(1);
	}
	int i,j,k;

	for(i = 0; i < NFILES; i++)
	for(j = 0; j < NCOOLINGSTEPS; j++)
	for(k = 0; k < NT; k++)
	{
		fscanf(input, "%f", &T23[i][j][k]);
		fscanf(input, "%f", &T24[i][j][k]);
		fscanf(input, "%f", &T34[i][j][k]);
	}
	fclose(input);
}

float average(float *arr, int L)
{
	int i;
	float av = 0;
	for(i = 0; i < L; i++)
	{
		av += arr[i];	
	}
	return av/L;
}

float error(float *arr, int L)
{
	int i;
	float ssq = 0;
	float av = average(arr, L);
	for(i = 0; i < L; i++)
	{
		ssq += arr[i]*arr[i];	 
	}	
	return sqrt(ssq/L - av*av);
}


void makeOutput(float Correlator[NCOOLINGSTEPS][NT], float Error[NCOOLINGSTEPS][NT])
{
	int i, dt;
	ofstream myfile;
	for(i = 0; i < NCOOLINGSTEPS; i++)
	{
		stringstream ss;
		ss << i;
		string name = "t24cool" + ss.str();
		myfile.open(name.c_str());
		
		for(dt = 0; dt < NT; dt++)
			myfile << dt << " " << Correlator[i][dt] << " " << Error[i][dt]  <<"\n";

		myfile.close();
	}
}
